Direct S3 Data Access

HTTPS vs s3

NASA Eartdata cloud provides two pathways for accessing data from the cloud. The first is via HTTPS. The other is through direct S3 bucket access.

Below are some considerations when selecting … Below are some advantages to choosing one or the other…

Objective

• Configure our notebook environment and retrieve temporary S3 credentials for in-region direct S3 bucket access
• Access a single HLS file
• Access and clip an HLS file to a region of interest
• Create an HLS time series data array

Let’s get started!

Import Required Packages

%matplotlib inline
import matplotlib.pyplot as plt
from datetime import datetime
import os
import subprocess
import requests
import boto3
import pandas as pd
import numpy as np
import xarray as xr
import rasterio as rio
from rasterio.session import AWSSession
from rasterio.plot import show
import rioxarray
import geopandas
import pyproj
from pyproj import Proj
from shapely.ops import transform
import geoviews as gv
from cartopy import crs
import hvplot.xarray
import holoviews as hv
gv.extension('bokeh', 'matplotlib')

Insert the credentials into our boto3 session and configure our rasterio environment for data access

Create a boto3 Session object using your temporary credentials. This Session is used to pass credentials and configuration to AWS so we can interact wit S3 objects from applicable buckets.

session = boto3.Session(aws_access_key_id=temp_creds_req['accessKeyId'], 
                        aws_secret_access_key=temp_creds_req['secretAccessKey'],
                        aws_session_token=temp_creds_req['sessionToken'],
                        region_name='us-west-2')

For this exercise, we are going to open up a context manager for the notebook using the rasterio.env module to store the required GDAL and AWS configurations we need to access the data in Earthdata Cloud. While the context manager is open (rio_env.__enter__()) we will be able to run the open or get data commands that would typically be executed within a with statement, thus allowing us to more freely interact with the data. We’ll close the context (rio_env.__exit__()) at the end of the notebook.

GDAL environment variables must be configured to access Earthdata Cloud data assets. Geospatial data access Python packages like rasterio and rioxarray depend on GDAL, leveraging GDAL’s “Virtual File Systems” to read remote files. GDAL has a lot of environment variables that control it’s behavior. Changing these settings can mean the difference being able to access a file or not. They can also have an impact on the performance.

rio_env = rio.Env(AWSSession(session),
                  GDAL_DISABLE_READDIR_ON_OPEN='EMPTY_DIR',
                  GDAL_HTTP_COOKIEFILE=os.path.expanduser('~/cookies.txt'),
                  GDAL_HTTP_COOKIEJAR=os.path.expanduser('~/cookies.txt'))
rio_env.__enter__()

List the available files in the data directory

for f in os.listdir('./data'):
    print(f)

S3_T12TGF_B04_Links.txt
nir_stack.vrt
S3_T12TGF_NIR_VSI_Links.txt
S3_T12TGF_Fmask_Links.txt
HTTPS_T12TGF_B8A_Links.txt
ne_w_agfields.geojson
red_stack.vrt
.ipynb_checkpoints
S3_T12TGF_RED_VSI_Links.txt
S3_T12TGF_B02_Links.txt
HTTPS_T12TGF_B05_Links.txt
S3_T12TGF_B05_Links.txt
fmask_stack.vrt
HTTPS_T12TGF_Fmask_Links.txt
S3_T12TGF_B8A_Links.txt
S3_T12TGF_FMASK_VSI_Links.txt
HTTPS_T12TGF_B04_Links.txt
HTTPS_T12TGF_B02_Links.txt

We will safe our list of links and a single link as Python objects for use later.

s3_links = open('./data/S3_T12TGF_B04_Links.txt').read().splitlines()
s3_links

s3_link = s3_links[0]
s3_link

Plot the HLS S3 object

hls_da.hvplot.image(x='x', y='y', cmap='fire', rasterize=True, width=800, height=600, colorbar=True)    # colormaps -> https://holoviews.org/user_guide/Colormaps.html

Unable to display output for mime type(s): 

We can print out the data value as a numpy array by typing .values

hls_da.values

Up to this point, we have not saved anything but metadata into memory. To save or load the data into memory we can call the .load() function.

hls_da_data = hls_da.load()
hls_da_data

<xarray.DataArray (y: 3660, x: 3660)>
array([[-9999, -9999, -9999, ...,  1527,  1440,  1412],
       [-9999, -9999, -9999, ...,  1493,  1476,  1407],
       [-9999, -9999, -9999, ...,  1466,  1438,  1359],
       ...,
       [-9999, -9999, -9999, ...,  1213,  1295,  1159],
       [-9999, -9999, -9999, ...,  1042,  1232,  1185],
       [-9999, -9999, -9999, ...,   954,  1127,  1133]], dtype=int16)
Coordinates:
  * x            (x) float64 7e+05 7e+05 7e+05 ... 8.097e+05 8.097e+05 8.097e+05
  * y            (y) float64 4.6e+06 4.6e+06 4.6e+06 ... 4.49e+06 4.49e+06
    spatial_ref  int64 0
Attributes:
    _FillValue:    -9999.0
    scale_factor:  0.0001
    add_offset:    0.0
    long_name:     Red

del(hls_da_data)

Read in our geojson file and transform its CRS

field = geopandas.read_file('./data/ne_w_agfields.geojson')

Let’s take a look at the bounding coordinate values.

field_shape = field.geometry[0]
field_shape.bounds

Note, the values above are in decimal degrees and represent the longitude and latitude for the lower left corner (-101.67271614074707, 41.04754380304359) and upper right corner (-101.65344715118408, 41.06213891056728) respectively.

Get the projection information from the HLS file

hls_proj = hls_da.rio.crs
hls_proj

Transform coordinates from lat lon (units = dd) to UTM (units = m)

geo_CRS = Proj('+proj=longlat +datum=WGS84 +no_defs', preserve_units=True)   # Source coordinate system of the ROI

project = pyproj.Transformer.from_proj(geo_CRS, hls_proj)                    # Set up the transformation

fsUTM = transform(project.transform, field_shape)
fsUTM.bounds

The coordinates for our feature have now been converted to UTM Zone 13N whether meters is the designated unit. Note the difference in the values between field_shape.bounds (in geographic) and fsUTM.bounds (in UTM projection).

Now we can clip our HLS file to our region of insterest!

Access and clip the HLS file

We can now use our transformed ROI bounding box to clip the HLS S3 object we accessed before. We’ll use the `rio.clip

hls_da_clip = rioxarray.open_rasterio(s3_link, chunks=True).squeeze('band', drop=True).rio.clip([fsUTM])

hls_da_clip

<xarray.DataArray (y: 56, x: 56)>
dask.array<astype, shape=(56, 56), dtype=int16, chunksize=(56, 56), chunktype=numpy.ndarray>
Coordinates:
  * y            (y) float64 4.551e+06 4.551e+06 ... 4.549e+06 4.549e+06
  * x            (x) float64 7.796e+05 7.796e+05 ... 7.812e+05 7.812e+05
    spatial_ref  int64 0
Attributes:
    scale_factor:  0.0001
    add_offset:    0.0
    long_name:     Red
    _FillValue:    -9999

hls_da_clip.hvplot.image(x = 'x', y = 'y', crs = 'EPSG:32613', cmap='fire', rasterize=True, width=800, height=600, colorbar=True)

Unable to display output for mime type(s): 

Read In and Process STAC Asset Links

In the previous section, we used the NASA CMR-STAC API to discover HLS assets the intersect with our search criteria, i.e., ROI, Date range, and collections. The search results were filtered and saved as text files by individual bands for each tile. We will read in the text files for tile T13TGF for the RED (L30: B04 & S30: B04), NIR (L30: B05 & S30: B8A), and Fmask bands.

[t for t in os.listdir('./data') if '.txt' in t]

red_s3_links = open('./data/S3_T12TGF_B04_Links.txt').read().splitlines()
red_s3_links

nir_bands = ['B05', 'B8A']
nir_link_text = [x for x in os.listdir('./data') if any(b in x for b in nir_bands) and 'S3' in x]
nir_s3_links = []
for file in nir_link_text:
    nir_s3_links.extend(open(f'./data/{file}').read().splitlines())
nir_s3_links

fmask_s3_links = open('./data/S3_T12TGF_Fmask_Links.txt').read().splitlines()
#fmask_s3_links

/vsis3/lp-prod-protected/HLSS30.015/HLS.S30.T13TGF.2020191T172901.v1.5.B04.tif

with open('./data/S3_T12TGF_RED_VSI_Links.txt', 'w') as f:
    links_vsi = [r.replace('s3://', '/vsis3/' ) + '\n' for r in red_s3_links]
    for link in links_vsi:
        f.write(link)

with open('./data/S3_T12TGF_NIR_VSI_Links.txt', 'w') as f:
    links_vsi = [r.replace('s3://', '/vsis3/' ) + '\n' for r in nir_s3_links]
    for link in links_vsi:
        f.write(link)

with open('./data/S3_T12TGF_FMASK_VSI_Links.txt', 'w') as f:
    links_vsi = [r.replace('s3://', '/vsis3/' ) + '\n' for r in fmask_s3_links]
    for link in links_vsi:
        f.write(link)

Read in geoJSON for subsetting

We will use the input geoJSON file to clip the source data to our desired region of interest.

field = geopandas.read_file('./data/ne_w_agfields.geojson')
fieldShape = field['geometry'][0]  

To clip the source data to our input feature boundary, we need to transform the feature boundary from its original WGS84 coordinate reference system to the projected reference system of the source HLS file (i.e., UTM Zone 13).

foa_url = red_s3_links[0]
with rio.open(foa_url) as src:
    hls_proj = src.crs.to_string()

hls_proj    

geo_CRS = Proj('+proj=longlat +datum=WGS84 +no_defs', preserve_units=True)   # Source coordinate system of the ROI
project = pyproj.Transformer.from_proj(geo_CRS, hls_proj)                    # Set up the transformation
fsUTM = transform(project.transform, fieldShape)

Direct S3 Data Access

#from dask.distributed import Client

#client = Client(n_workers=2)
#client

build_red_vrt = f"gdalbuildvrt ./data/red_stack.vrt -separate -input_file_list ./data/S3_T12TGF_RED_VSI_Links.txt --config AWS_ACCESS_KEY_ID {temp_creds_req['accessKeyId']} --config AWS_SECRET_ACCESS_KEY {temp_creds_req['secretAccessKey']} --config AWS_SESSION_TOKEN {temp_creds_req['sessionToken']} --config GDAL_DISABLE_READDIR_ON_OPEN TRUE"
#build_red_vrt    # !!! BEWARE, removing the # on this line will print your temporary S3 credentials.

%%time

subprocess.call(build_red_vrt, shell=True)

build_nir_vrt = f"gdalbuildvrt ./data/nir_stack.vrt -separate -input_file_list ./data/S3_T12TGF_NIR_VSI_Links.txt --config AWS_ACCESS_KEY_ID {temp_creds_req['accessKeyId']} --config AWS_SECRET_ACCESS_KEY {temp_creds_req['secretAccessKey']} --config AWS_SESSION_TOKEN {temp_creds_req['sessionToken']} --config GDAL_DISABLE_READDIR_ON_OPEN TRUE"
subprocess.call(build_nir_vrt, shell=True)

build_fmask_vrt = f"gdalbuildvrt ./data/fmask_stack.vrt -separate -input_file_list ./data/S3_T12TGF_FMASK_VSI_Links.txt --config AWS_ACCESS_KEY_ID {temp_creds_req['accessKeyId']} --config AWS_SECRET_ACCESS_KEY {temp_creds_req['secretAccessKey']} --config AWS_SESSION_TOKEN {temp_creds_req['sessionToken']} --config GDAL_DISABLE_READDIR_ON_OPEN TRUE"
subprocess.call(build_fmask_vrt, shell=True)

Reading in an HLS time series

We can now read the VRT files into our Python session. A drawback of reading VRTs into Python is that the time coordinate variable needs to be contructed. Below we not only read in the VRT file using rioxarray, but we also repurpose the band variable, which is generated automatically, to hold out time information.

%%time

chunks=dict(band=1, x=1024, y=1024)
#chunks=dict(band=1, x=512, y=512)
red = rioxarray.open_rasterio('./data/red_stack.vrt', chunks=chunks)                    # Read in VRT
red = red.rename({'band':'time'})                                                       # Rename the 'band' coordinate variable to 'time' 
red['time'] = [datetime.strptime(x.split('.')[-5], '%Y%jT%H%M%S') for x in links_vsi]   # Extract the time information from the input file names and assign them to the time coordinate variable
red = red.sortby('time')                                                                # Sort by the time coordinate variable
red

red.time

red_clip = red.rio.clip([fsUTM]).persist()
red_clip

red_clip.hvplot.image(x='x', y='y', width=800, height=600, colorbar=True, cmap='Reds').opts(clim=(0.0, red_clip.values.max()))

Read in the NIR and Fmask VRT files

%%time
chunks=dict(band=1, x=1024, y=1024)
nir = rioxarray.open_rasterio('./data/nir_stack.vrt', chunks=chunks)                    # Read in VRT
nir = nir.rename({'band':'time'})                                                       # Rename the 'band' coordinate variable to 'time' 
nir['time'] = [datetime.strptime(x.split('.')[-5], '%Y%jT%H%M%S') for x in links_vsi]   # Extract the time information from the input file names and assign them to the time coordinate variable
nir = nir.sortby('time')                                                                # Sort by the time coordinate variable
nir

%%time
chunks=dict(band=1, x=1024, y=1024)
fmask = rioxarray.open_rasterio('./data/fmask_stack.vrt', chunks=chunks)                    # Read in VRT
fmask = fmask.rename({'band':'time'})                                                       # Rename the 'band' coordinate variable to 'time' 
fmask['time'] = [datetime.strptime(x.split('.')[-5], '%Y%jT%H%M%S') for x in links_vsi]     # Extract the time information from the input file names and assign them to the time coordinate variable
fmask = fmask.sortby('time')                                                                # Sort by the time coordinate variable
fmask

Create an xarray dataset

We will now combine the RED, NIR, and Fmask arrays into a dataset and create/add a new NDVI variable.

hls_ndvi = xr.Dataset({'red': red, 'nir': nir, 'fmask': fmask, 'ndvi': (nir - red) / (nir + red)})
hls_ndvi

Above, we created a new NDVI variable. Now, we will clip and plot our results.

ndvi_clip = hls_ndvi.ndvi.rio.clip([fsUTM]).persist()
ndvi_clip

ndvi_clip.hvplot.image(x='x', y='y', groupby='time', width=800, height=600, colorbar=True, cmap='YlGn').opts(clim=(0.0, 1.0))

Apply quality filter

We want to keep NDVI data values where Fmask equals 0 (no clouds, no cloud shadow, no snow/ice, no water.

ndvi_clip_filter = hls_ndvi.ndvi.where(fmask==0, np.nan).rio.clip([fsUTM]).persist()

ndvi_clip_filter.hvplot.image(x='x', y='y', groupby='time', width=800, height=600, colorbar=True, cmap='YlGn').opts(clim=(0.0, 1.0))

ndvi_clip_filter.groupby('time.month').mean('time').hvplot.image(x = 'x', y = 'y', crs = hls_proj, groupby='month', cmap='YlGn', width=800, height=600, colorbar=True).opts(clim=(0.0, 1.0))

rio_env.__exit__()